SSVEP reward - EEG results
Ivan Grahek & Antonio Schettino
2019-07-29
About the code
Experiment: FSAReward (Ivan Grahek, Antonio Schettino, Gilles Pourtois, Ernst Koster, & Søren Andersen) (*: co-first authors) Code written by: Ivan Grahek & Antonio Schettino (2016-2018) Description: Code for the analysis of EEG data for Experiment 1 of the SSVEP - reward project.
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
EEG
Topography & spectra
Import data
# # Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# clear the environment
rm(list=ls())
# clear the console
cat("\014")
#load packages and install them if they're not installed
if (!require("pacman")) install.packages("pacman")
pacman::p_load(reshape2,yarrr,BayesFactor,plyr,ez,schoRsch,brms,knitr,broom, brmstools, BEST, tidyverse, here, HDInterval)
# set seed
set.seed(42)
# Set working directory
#setwd(here())
# import data
data.raw = read.csv(file = here("EEG_preprocessing/movement","grandAverage_amplitudes.csv"),header=TRUE,na.strings="NaN")
# Rename columns
colnames(data.raw) = c("Subject","Frequency","BslnRedAttendedNomov","BslnBlueAttendedNomov","AcqRedAttendedNomov","AcqBlueAttendedNomov","ExtRedAttendedNomov","ExtBlueAttendedNomov","BslnRedAttendedMov","BslnBlueAttendedMov","AcqRedAttendedMov","AcqBlueAttendedMov","ExtRedAttendedMov","ExtBlueAttendedMov")
# # Reshape to long format - take only no movement trials
data = melt(data.raw,id.vars=c("Subject","Frequency"),
measure.vars=c("BslnRedAttendedNomov","BslnBlueAttendedNomov","AcqRedAttendedNomov","AcqBlueAttendedNomov","ExtRedAttendedNomov","ExtBlueAttendedNomov","BslnRedAttendedMov","BslnBlueAttendedMov","AcqRedAttendedMov","AcqBlueAttendedMov","ExtRedAttendedMov","ExtBlueAttendedMov"),
variable.name="Condition",value.name="Amplitude")
# Sort the new dataframe by participant name
data = data[order(data$Subject),]
# Split the variable Condition based on capital letters
data$Condition = gsub("(?!^)(?=[[:upper:]])", " ", data$Condition, perl=T)
# Split the variable condition into multiple variables
Conditions = colsplit(data$Condition, pattern="\\s+",names = c('ExpPhase', 'ColorMoved',"attended","no","moved"))
# Add the variable defining which color is rewarded based on the participant number
data$RewardedColor = ifelse(data$Subject%%2==0,"Blue","Red") # if participant number is even, blue was rewarded
# Add the Conditions needed to the dataset
data$ExpPhase = Conditions[,1]
data$AttendedColor = Conditions[,2]
data$Movement = Conditions[,4]
# Switch the Frequency to the color
data$RecordedFrequency = ifelse(data$Frequency==10,"Blue","Red") # if the recorded frequency is 10Hz assign Blue (color flickering at 10Hz), otherwise assign Red (color flickering at 12Hz)
# Make a new condition based on the attended color and the rewarded color
data$Condition = ifelse(data$AttendedColor==data$RewardedColor, "High_Rew","Low_Rew")
# Make a new condition based on the attended color and the recorded frequency
data$Attention = ifelse(data$AttendedColor==data$RecordedFrequency, "Att","NotAtt")
# Make a new condition based the Condition and the Attention
data$RecordingAndCondition = with(data, paste0(Condition,"_",Attention))
# Select variables which we want to keep
data = subset(data, select=c("Subject","RewardedColor","ExpPhase","AttendedColor","Condition","RecordedFrequency","Attention","RecordingAndCondition","Movement","Amplitude"))
# Sort the data
data = data[with(data, order(Subject)), ]
# Make a new variable with mean amplitude across all conditions for each participant and each frequency (normalization)
data = ddply(data,.(Subject,RecordedFrequency),transform,
MeanAmplitude = mean(Amplitude[ExpPhase=="Bsln"],na.rm=TRUE),
SDAmplitude = sd(Amplitude,na.rm=TRUE))
# Divide amplitudes in each Subject, Frequency, and Condition by the Mean Amplitude
data$Amplitude = data$Amplitude/data$MeanAmplitude
# Calculate the attention indexes - Selectivity (attended-unattended) & total enhancement (attended+unattended) (Andersen & Muller, 2010, PNAS)
data.diff = ddply(data, .(Subject,ExpPhase,Condition), transform, Selectivity = Amplitude[Attention=="Att"]-Amplitude[Attention=="NotAtt"],TotalEnhancement=Amplitude[Attention=="Att"]+Amplitude[Attention=="NotAtt"])
# Delete the Attention column and rows which are not necessary (indexes repeated twice)
data.diff = subset(data.diff,Attention=="Att") #keep only Att as it is equal to NotAtt
data.diff$Attention = NULL #drop the Attention column
# Sort the data
data.diff$ExpPhase = factor(data.diff$ExpPhase, levels = c("Bsln","Acq","Ext"))
data.diff = data.diff[order(data.diff$Subject,data.diff$Condition,data.diff$ExpPhase),]
# Calculate the reward index - High reward minus Low reward
data.reward = ddply(data, .(Subject,ExpPhase,Attention), transform, Reward = Amplitude[Condition=="High_Rew"]-Amplitude[Condition=="Low_Rew"])
# Delete the Attention column and rows which are not necessary (indexes repeated twice)
data.reward = subset(data.reward,Condition=="High_Rew") #keep only Att as it is equal to NotAtt
data.reward$Condition = NULL #drop the Condition column
# Sort the data
data.reward$ExpPhase = factor(data.reward$ExpPhase, levels = c("Bsln","Acq","Ext"))
data.reward = data.reward[order(data.reward$Subject,data.reward$Attention,data.reward$ExpPhase),]
# Convert variables to be used in analyses into factors
data[c("Subject", "Condition","ExpPhase", "RewardedColor", "Attention", "RecordingAndCondition")] =
lapply(data[c("Subject", "Condition","ExpPhase", "RewardedColor", "Attention", "RecordingAndCondition")], factor)
data.diff[c("Subject", "Condition","ExpPhase", "RecordingAndCondition")] =
lapply(data.diff[c("Subject", "Condition","ExpPhase", "RecordingAndCondition")], factor)
# Save the final data into a new variable
data_all = dataAll trials
Means - raw data
summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))
names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")
summary$Attention = recode(summary$Attention,
"Att" = "Attended",
"NotAtt" = "Unattended")
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL
kable(summary, caption = "Amplitudes per condition")| Attention | Reward phase | Reward probability | Amplitude |
|---|---|---|---|
| Attended | Baseline | High | 1.13 [ 0.92 1.52 ] |
| Attended | Baseline | Low | 1.13 [ 0.86 1.52 ] |
| Attended | Acquisition | High | 1.16 [ 0.8 1.6 ] |
| Attended | Acquisition | Low | 1.13 [ 0.76 1.71 ] |
| Attended | Extinction | High | 1.13 [ 0.61 1.61 ] |
| Attended | Extinction | Low | 1.13 [ 0.59 1.84 ] |
| Unattended | Baseline | High | 0.87 [ 0.47 1.17 ] |
| Unattended | Baseline | Low | 0.87 [ 0.49 1.11 ] |
| Unattended | Acquisition | High | 0.91 [ 0.54 1.38 ] |
| Unattended | Acquisition | Low | 0.89 [ 0.5 1.28 ] |
| Unattended | Extinction | High | 0.88 [ 0.48 1.23 ] |
| Unattended | Extinction | Low | 0.91 [ 0.44 1.42 ] |
Plots - raw data
# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"
# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]
plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.2,2.2), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
}
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.alltrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.alltrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.alltrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.alltrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.alltrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.alltrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.alltrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.alltrials")
bR2.full.EEG = readRDS("bR2.full.EEG.alltrials")Model comparison with WAIC
print(compare.threefactors.EEG.waic)## WAIC SE
## null 65.34 62.67
## expphase 16.17 59.31
## attention -341.85 70.91
## phaseANDattention -428.76 67.19
## phaseANDattention_interaction -424.98 67.36
## rewardTimesPhasePlusAtt -647.37 74.90
## full -638.31 74.97
## null - expphase 49.17 16.02
## null - attention 407.18 38.39
## null - phaseANDattention 494.10 41.45
## null - phaseANDattention_interaction 490.31 41.38
## null - rewardTimesPhasePlusAtt 712.71 49.62
## null - full 703.65 49.28
## expphase - attention 358.01 42.85
## expphase - phaseANDattention 444.93 38.51
## expphase - phaseANDattention_interaction 441.14 38.39
## expphase - rewardTimesPhasePlusAtt 663.54 48.36
## expphase - full 654.48 47.97
## attention - phaseANDattention 86.92 24.12
## attention - phaseANDattention_interaction 83.13 24.23
## attention - rewardTimesPhasePlusAtt 305.53 39.74
## attention - full 296.46 39.74
## phaseANDattention - phaseANDattention_interaction -3.79 1.34
## phaseANDattention - rewardTimesPhasePlusAtt 218.61 35.06
## phaseANDattention - full 209.55 35.07
## phaseANDattention_interaction - rewardTimesPhasePlusAtt 222.40 35.06
## phaseANDattention_interaction - full 213.33 35.01
## rewardTimesPhasePlusAtt - full -9.06 3.03Bayesian R squared
Null model
print(bR2.null.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.07786497 0.01721677 0.04563591 0.1124821Experiment phase
print(bR2.expphase.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.1424397 0.02004885 0.1036544 0.1816891Attention model
print(bR2.attention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.3975517 0.01994732 0.3571852 0.4352368Phase and attention model
print(bR2.phaseANDattention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4637211 0.018479 0.4260748 0.4983349Phase and attention interaction model
print(bR2.phaseANDattention_interaction.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4640534 0.0185034 0.4263085 0.4983363Reward times phase plus attention model
print(bR2.rewardTimesPhasePlusAtt.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.5882896 0.01484724 0.5579976 0.6158518Full model
print(bR2.full.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.5881777 0.01486025 0.557527 0.6157646Checking the best model
Plotting the chains
# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 1.125 0.02 1.091 1.16
## 2 ConditionLow_Rew -0.001 0.03 -0.054 0.05
## 3 ExpPhaseAcq 0.038 0.03 -0.003 0.08
## 4 ExpPhaseExt 0.006 0.03 -0.039 0.05
## 5 AttentionNotAtt -0.249 0.02 -0.289 -0.21
## 6 ConditionLow_Rew:ExpPhaseAcq -0.028 0.02 -0.068 0.01
## 7 ConditionLow_Rew:ExpPhaseExt 0.015 0.02 -0.026 0.06Plotting the best model
post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")
# Calculate posteriors for each condition
################################################ Baseline ####
##################### Attended
######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
##################### Not Attended
######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]]
################################################ Acquistion
##################### Attended
######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]]
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
##################### Not Attended
######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
################################################ Extinction
##################### Attended
######### High reward
Extinction_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]]
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]]
##################### Not Attended
######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]] # Data from the model
# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")
posterior_conditions$Attention = recode(posterior_conditions$Attention,
"Attended" = "Target",
"NotAttended" = "Distractor")
posterior_conditions$Attention = as.factor(posterior_conditions$Attention )
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward probability",
"Low" = "Low reward probability")
posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)
names(posterior_conditions)[4] = "Amplitude"
# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Unattended Baseline High Reward",
"Attended Baseline Low Reward",
"Unattended Baseline Low Reward",
"Attended Acquisition High Reward",
"Unattended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Unattended Acquisition Low Reward",
"Attended Extinction High Reward",
"Unattended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_High_Attended)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_High_Attended)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward probability",
"Low Reward" = "Low reward probability")
posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)
posterior_means$Attention = recode(posterior_means$Attention,
"Attended" = "Target",
"Unattended" = "Distractor")
posterior_means$Attention = as.factor(posterior_means$Attention)
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)
# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion
dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"
# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
"Acq" = "Training",
"Bsln" = "Baseline",
"Ext" = "Test")
dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
"High_Rew" = "High reward probability",
"Low_Rew" = "Low reward probability")
dataPlot$Attention = recode(dataPlot$Attention,
"Att" = "Target",
"NotAtt" = "Distractor")
##### Plotting targets and distractors separately #####
# Attended
dataAmplitudePlot=subset(dataPlot,Attention=="Target")
posterior_means_plot=subset(posterior_means,Attention=="Target")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Target")
# tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Unattended
dataAmplitudePlot=subset(dataPlot,Attention=="Distractor")
posterior_means_plot=subset(posterior_means,Attention=="Distractor")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Distractor")
# tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
##### Plotting high and low separately #####
# High
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="High reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Low
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="Low reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
####### Plotting everything together #####
# Order for plotting
# dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor"))
# posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor"))
# posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor"))
#
# dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)
ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(`Reward Probability`~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
theme(
axis.line = element_line(size=2, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 24, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 22),
axis.text.x=element_text(colour="black", size = 22),
axis.text.y=element_text(colour="black", size = 22),
axis.title=element_text(size=24,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 24))
# dev.off()Table of means across conditions
# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Attended Baseline Low Reward",
"Attended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Attended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Baseline High Reward",
"Unattended Baseline Low Reward",
"Unattended Acquisition High Reward",
"Unattended Acquisition Low Reward",
"Unattended Extinction High Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Attention | Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|---|
| Attended | Baseline | High Reward | 1.12 [ 1.08 1.17 ] |
| Attended | Baseline | Low Reward | 1.12 [ 1.08 1.17 ] |
| Attended | Acquisition | High Reward | 1.16 [ 1.11 1.22 ] |
| Attended | Acquisition | Low Reward | 1.13 [ 1.07 1.19 ] |
| Attended | Extinction | High Reward | 1.13 [ 1.07 1.19 ] |
| Attended | Extinction | Low Reward | 0.9 [ 0.84 0.96 ] |
| Unattended | Baseline | High Reward | 0.88 [ 0.83 0.92 ] |
| Unattended | Baseline | Low Reward | 0.87 [ 0.83 0.92 ] |
| Unattended | Acquisition | High Reward | 0.91 [ 0.86 0.97 ] |
| Unattended | Acquisition | Low Reward | 0.89 [ 0.83 0.94 ] |
| Unattended | Extinction | High Reward | 0.88 [ 0.82 0.94 ] |
| Unattended | Extinction | Low Reward | 0.9 [ 0.84 0.96 ] |
| ### Inference | about the best | model |
Attended vs. unattended
Check the difference between attended and not attended in baseline high rewarded
Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in baseline low rewarded
Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition high rewarded
Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition low rewarded
Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction high rewarded
Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction low rewarded
Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Comparison between phases
Check the difference between baseline and acquisition in high reward attended
Diff_Bsln_Acq_High_Att = Acquisition_High_Attended - Baseline_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward attended
Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward not attended
Diff_Bsln_Acq_High_NotAtt = Acquisition_High_NotAttended - Baseline_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward not attended
Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward attended
Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward attended
Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward not attended
Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward not attended
Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Baseline difference
Check the difference between high and low reward in baseline attended
Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between high and low reward in baseline not attended
Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward unattended vs. attended
Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")## [1] "Mean = 0 [ 0 0 ]"Only movement trials
# Take only the movement trials
data = subset(data_all, Movement=="Mov")Means - raw data
summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))
names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")
summary$Attention = recode(summary$Attention,
"Att" = "Attended",
"NotAtt" = "Unattended")
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL
kable(summary, caption = "Amplitudes per condition")| Attention | Reward phase | Reward probability | Amplitude |
|---|---|---|---|
| Attended | Baseline | High | 1.1 [ 0.77 1.39 ] |
| Attended | Baseline | Low | 1.12 [ 0.86 1.52 ] |
| Attended | Acquisition | High | 1.17 [ 0.83 1.71 ] |
| Attended | Acquisition | Low | 1.12 [ 0.8 1.48 ] |
| Attended | Extinction | High | 1.15 [ 0.61 1.6 ] |
| Attended | Extinction | Low | 1.11 [ 0.59 1.84 ] |
| Unattended | Baseline | High | 0.87 [ 0.5 1.17 ] |
| Unattended | Baseline | Low | 0.87 [ 0.66 1.07 ] |
| Unattended | Acquisition | High | 0.9 [ 0.59 1.21 ] |
| Unattended | Acquisition | Low | 0.87 [ 0.69 1.15 ] |
| Unattended | Extinction | High | 0.86 [ 0.44 1.1 ] |
| Unattended | Extinction | Low | 0.9 [ 0.44 1.28 ] |
Plots - raw data
# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"
# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]
plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.2,2.2), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
}
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.movementtrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.movementtrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.movementtrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.movementtrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.movementtrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.movementtrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.movementtrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.movementtrials")
bR2.full.EEG = readRDS("bR2.full.EEG.movementtrials")Model comparison with WAIC
print(compare.threefactors.EEG.waic)## WAIC SE
## null 24.68 48.40
## expphase 17.95 46.29
## attention -158.66 55.02
## phaseANDattention -181.20 52.08
## phaseANDattention_interaction -177.80 51.70
## rewardTimesPhasePlusAtt -290.73 58.73
## full -282.50 58.65
## null - expphase 6.73 8.46
## null - attention 183.35 27.25
## null - phaseANDattention 205.88 27.81
## null - phaseANDattention_interaction 202.48 27.79
## null - rewardTimesPhasePlusAtt 315.41 31.89
## null - full 307.18 31.39
## expphase - attention 176.61 29.33
## expphase - phaseANDattention 199.15 26.68
## expphase - phaseANDattention_interaction 195.75 26.66
## expphase - rewardTimesPhasePlusAtt 308.68 32.33
## expphase - full 300.45 31.80
## attention - phaseANDattention 22.54 14.45
## attention - phaseANDattention_interaction 19.13 14.45
## attention - rewardTimesPhasePlusAtt 132.06 26.06
## attention - full 123.84 26.10
## phaseANDattention - phaseANDattention_interaction -3.40 1.04
## phaseANDattention - rewardTimesPhasePlusAtt 109.53 25.57
## phaseANDattention - full 101.30 25.46
## phaseANDattention_interaction - rewardTimesPhasePlusAtt 112.93 25.75
## phaseANDattention_interaction - full 104.70 25.61
## rewardTimesPhasePlusAtt - full -8.23 4.32Bayesian R squared
Null model
print(bR2.null.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.07035095 0.02712311 0.01932853 0.125171Experiment phase
print(bR2.expphase.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.13111 0.03150338 0.07034581 0.1937305Attention model
print(bR2.attention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.3792606 0.03124579 0.3159362 0.4373363Phase and attention model
print(bR2.phaseANDattention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4455475 0.03195549 0.3795457 0.5038668Phase and attention interaction model
print(bR2.phaseANDattention_interaction.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4455208 0.03192962 0.3799893 0.5043975Reward times phase plus attention model
print(bR2.rewardTimesPhasePlusAtt.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.6015452 0.02468267 0.5491233 0.6460264Full model
print(bR2.full.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.6020692 0.02463371 0.5501192 0.6467771Checking the best model
Plotting the chains
# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 1.113 0.03 1e+00 1.156
## 2 ConditionLow_Rew 0.007 0.04 -5e-02 0.065
## 3 ExpPhaseAcq 0.047 0.03 2e-04 0.094
## 4 ExpPhaseExt 0.018 0.03 -4e-02 0.071
## 5 AttentionNotAtt -0.249 0.02 -3e-01 -0.211
## 6 ConditionLow_Rew:ExpPhaseAcq -0.047 0.03 -1e-01 0.009
## 7 ConditionLow_Rew:ExpPhaseExt -0.007 0.03 -6e-02 0.049Plotting the best model
post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")
# Calculate posteriors for each condition
################################################ Baseline ####
##################### Attended
######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
##################### Not Attended
######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]]
################################################ Acquistion
##################### Attended
######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]]
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
##################### Not Attended
######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
################################################ Extinction
##################### Attended
######### High reward
Extinction_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]]
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]]
##################### Not Attended
######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]] # make a data frame
posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")
posterior_conditions$Attention = recode(posterior_conditions$Attention,
"Attended" = "Attended",
"NotAttended" = "Unattended")
names(posterior_conditions)[4] = "Amplitude"
#order
#dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
#dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]
plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Attended")}
if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Unattended")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Reward Phase` + `Reward Probability`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.7,1.3), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
sortx = "sequential",
bty="l", # plot box type
back.col="white") # background, color
}
plottingConditions = c("High","Low" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="High"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="High")}
if(plottingConditions[i]=="Low"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="Low")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Attention` + `Reward Phase`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.7,1.3), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
sortx = "sequential",
bty="l", # plot box type
back.col="white") # background, color
}
# Data from the model
# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")
posterior_conditions$Attention = recode(posterior_conditions$Attention,
"Attended" = "Target",
"NotAttended" = "Distractor")
posterior_conditions$Attention = as.factor(posterior_conditions$Attention )
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward probability",
"Low" = "Low reward probability")
posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)
names(posterior_conditions)[4] = "Amplitude"
# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Unattended Baseline High Reward",
"Attended Baseline Low Reward",
"Unattended Baseline Low Reward",
"Attended Acquisition High Reward",
"Unattended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Unattended Acquisition Low Reward",
"Attended Extinction High Reward",
"Unattended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_High_Attended)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_High_Attended)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward probability",
"Low Reward" = "Low reward probability")
posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)
posterior_means$Attention = recode(posterior_means$Attention,
"Attended" = "Target",
"Unattended" = "Distractor")
posterior_means$Attention = as.factor(posterior_means$Attention)
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)
# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion
dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"
# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
"Acq" = "Training",
"Bsln" = "Baseline",
"Ext" = "Test")
dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
"High_Rew" = "High reward probability",
"Low_Rew" = "Low reward probability")
dataPlot$Attention = recode(dataPlot$Attention,
"Att" = "Target",
"NotAtt" = "Distractor")
##### Plotting targets and distractors separately #####
# Attended
dataAmplitudePlot=subset(dataPlot,Attention=="Target")
posterior_means_plot=subset(posterior_means,Attention=="Target")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Target")
# tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Unattended
dataAmplitudePlot=subset(dataPlot,Attention=="Distractor")
posterior_means_plot=subset(posterior_means,Attention=="Distractor")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Distractor")
# tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
##### Plotting high and low separately #####
# High
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="High reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Low
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="Low reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
####### Plotting everything together #####
# Order for plotting
# dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor"))
# posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor"))
# posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor"))
#
# dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)
ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(`Reward Probability`~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
theme(
axis.line = element_line(size=2, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 24, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 22),
axis.text.x=element_text(colour="black", size = 22),
axis.text.y=element_text(colour="black", size = 22),
axis.title=element_text(size=24,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 24))
# dev.off()Table of means across conditions
# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Attended Baseline Low Reward",
"Attended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Attended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Baseline High Reward",
"Unattended Baseline Low Reward",
"Unattended Acquisition High Reward",
"Unattended Acquisition Low Reward",
"Unattended Extinction High Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Attention | Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|---|
| Attended | Baseline | High Reward | 1.11 [ 1.06 1.17 ] |
| Attended | Baseline | Low Reward | 1.12 [ 1.07 1.17 ] |
| Attended | Acquisition | High Reward | 1.16 [ 1.09 1.22 ] |
| Attended | Acquisition | Low Reward | 1.12 [ 1.06 1.18 ] |
| Attended | Extinction | High Reward | 1.13 [ 1.06 1.2 ] |
| Attended | Extinction | Low Reward | 0.88 [ 0.82 0.95 ] |
| Unattended | Baseline | High Reward | 0.86 [ 0.81 0.91 ] |
| Unattended | Baseline | Low Reward | 0.87 [ 0.83 0.92 ] |
| Unattended | Acquisition | High Reward | 0.91 [ 0.85 0.97 ] |
| Unattended | Acquisition | Low Reward | 0.87 [ 0.82 0.92 ] |
| Unattended | Extinction | High Reward | 0.88 [ 0.82 0.95 ] |
| Unattended | Extinction | Low Reward | 0.88 [ 0.82 0.95 ] |
| ### Inference | about the best | model |
Attended vs. unattended
Check the difference between attended and not attended in baseline high rewarded
Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in baseline low rewarded
Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition high rewarded
Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition low rewarded
Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction high rewarded
Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction low rewarded
Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Comparison between phases
Check the difference between baseline and acquisition in high reward attended
Diff_Bsln_Acq_High_Att = Baseline_High_Attended - Acquisition_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward attended
Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward not attended
Diff_Bsln_Acq_High_NotAtt = Baseline_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward not attended
Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward attended
Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward attended
Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward not attended
Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward not attended
Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Baseline difference
Check the difference between high and low reward in baseline attended
Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between high and low reward in baseline not attended
Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward unattended vs. attended
Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")## [1] "Mean = 0 [ 0 0 ]"Only no movement trials
# Take only the movement trials
data = subset(data_all, Movement=="Nomov")Means - raw data
summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))
names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")
summary$Attention = recode(summary$Attention,
"Att" = "Attended",
"NotAtt" = "Unattended")
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL
kable(summary, caption = "Amplitudes per condition")| Attention | Reward phase | Reward probability | Amplitude |
|---|---|---|---|
| Attended | Baseline | High | 1.16 [ 0.99 1.52 ] |
| Attended | Baseline | Low | 1.13 [ 0.9 1.44 ] |
| Attended | Acquisition | High | 1.16 [ 0.8 1.52 ] |
| Attended | Acquisition | Low | 1.15 [ 0.76 1.71 ] |
| Attended | Extinction | High | 1.12 [ 0.64 1.61 ] |
| Attended | Extinction | Low | 1.15 [ 0.85 1.97 ] |
| Unattended | Baseline | High | 0.87 [ 0.62 1.22 ] |
| Unattended | Baseline | Low | 0.87 [ 0.6 1.27 ] |
| Unattended | Acquisition | High | 0.93 [ 0.54 1.38 ] |
| Unattended | Acquisition | Low | 0.9 [ 0.62 1.43 ] |
| Unattended | Extinction | High | 0.89 [ 0.5 1.23 ] |
| Unattended | Extinction | Low | 0.92 [ 0.5 1.42 ] |
Plots - raw data
# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"
# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]
plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.2,2.2), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
}
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.nomovementtrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.nomovementtrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.nomovementtrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.nomovementtrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.nomovementtrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.nomovementtrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.nomovementtrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.nomovementtrials")
bR2.full.EEG = readRDS("bR2.full.EEG.nomovementtrials")Model comparison with WAIC
print(compare.threefactors.EEG.waic)## WAIC SE
## null 63.46 46.48
## expphase 56.87 44.87
## attention -122.28 52.10
## phaseANDattention -146.66 49.31
## phaseANDattention_interaction -143.44 49.58
## rewardTimesPhasePlusAtt -218.42 54.95
## full -208.38 55.07
## null - expphase 6.59 9.33
## null - attention 185.74 25.05
## null - phaseANDattention 210.12 27.03
## null - phaseANDattention_interaction 206.90 26.65
## null - rewardTimesPhasePlusAtt 281.88 32.63
## null - full 271.85 32.05
## expphase - attention 179.15 28.01
## expphase - phaseANDattention 203.53 24.98
## expphase - phaseANDattention_interaction 200.31 24.62
## expphase - rewardTimesPhasePlusAtt 275.29 32.12
## expphase - full 265.26 31.59
## attention - phaseANDattention 24.38 16.42
## attention - phaseANDattention_interaction 21.16 16.17
## attention - rewardTimesPhasePlusAtt 96.14 23.60
## attention - full 86.11 23.17
## phaseANDattention - phaseANDattention_interaction -3.22 2.17
## phaseANDattention - rewardTimesPhasePlusAtt 71.76 21.80
## phaseANDattention - full 61.72 21.70
## phaseANDattention_interaction - rewardTimesPhasePlusAtt 74.98 21.86
## phaseANDattention_interaction - full 64.94 21.54
## rewardTimesPhasePlusAtt - full -10.03 2.60Bayesian R squared
Null model
print(bR2.null.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.08451262 0.02800565 0.03170381 0.1417311Experiment phase
print(bR2.expphase.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.1375776 0.0313894 0.07685781 0.1993291Attention model
print(bR2.attention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4010355 0.03056095 0.338263 0.457752Phase and attention model
print(bR2.phaseANDattention.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4569734 0.02908633 0.3965695 0.5105794Phase and attention interaction model
print(bR2.phaseANDattention_interaction.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4586185 0.02890073 0.3986922 0.5119207Reward times phase plus attention model
print(bR2.rewardTimesPhasePlusAtt.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.5695265 0.02575516 0.5153769 0.615958Full model
print(bR2.full.EEG)## Estimate Est.Error Q2.5 Q97.5
## R2 0.5690858 0.02584686 0.5148487 0.6155357Checking the best model
Plotting the chains
# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 1.138 0.02 1.10 1.18
## 2 ConditionLow_Rew -0.012 0.04 -0.07 0.05
## 3 ExpPhaseAcq 0.029 0.03 -0.02 0.08
## 4 ExpPhaseExt -0.007 0.03 -0.06 0.05
## 5 AttentionNotAtt -0.247 0.03 -0.29 -0.20
## 6 ConditionLow_Rew:ExpPhaseAcq -0.007 0.04 -0.07 0.05
## 7 ConditionLow_Rew:ExpPhaseExt 0.039 0.04 -0.02 0.10Plotting the best model
post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")
# Calculate posteriors for each condition
################################################ Baseline ####
##################### Attended
######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
##################### Not Attended
######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]]
################################################ Acquistion
##################### Attended
######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]]
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
##################### Not Attended
######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseAcq"]]
################################################ Extinction
##################### Attended
######### High reward
Extinction_High_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]]
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]]
##################### Not Attended
######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_AttentionNotAtt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ConditionLow_Rew:ExpPhaseExt"]] # make a data frame
posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")
posterior_conditions$Attention = recode(posterior_conditions$Attention,
"Attended" = "Attended",
"NotAttended" = "Unattended")
names(posterior_conditions)[4] = "Amplitude"
#order
#dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
#dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]
plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Attended")}
if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Unattended")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Reward Phase` + `Reward Probability`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.7,1.3), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
sortx = "sequential",
bty="l", # plot box type
back.col="white") # background, color
}
plottingConditions = c("High","Low" )
for (i in 1:length(plottingConditions)){
if(plottingConditions[i]=="High"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="High")}
if(plottingConditions[i]=="Low"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="Low")}
# Pirate plot
pirateplot(formula = Amplitude ~ `Attention` + `Reward Phase`, # dependent~independent variables
data=dataAmplitudePlot, # data frame
main=plottingConditions[i], # main title
ylim=c(0.7,1.3), # y-axis: limits
ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
theme=0, # preset theme (0: use your own)
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
sortx = "sequential",
bty="l", # plot box type
back.col="white") # background, color
}
# Data from the model
# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")
posterior_conditions$Attention = recode(posterior_conditions$Attention,
"Attended" = "Target",
"NotAttended" = "Distractor")
posterior_conditions$Attention = as.factor(posterior_conditions$Attention )
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward probability",
"Low" = "Low reward probability")
posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)
names(posterior_conditions)[4] = "Amplitude"
# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Unattended Baseline High Reward",
"Attended Baseline Low Reward",
"Unattended Baseline Low Reward",
"Attended Acquisition High Reward",
"Unattended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Unattended Acquisition Low Reward",
"Attended Extinction High Reward",
"Unattended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_High_Attended)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_High_Attended)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward probability",
"Low Reward" = "Low reward probability")
posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)
posterior_means$Attention = recode(posterior_means$Attention,
"Attended" = "Target",
"Unattended" = "Distractor")
posterior_means$Attention = as.factor(posterior_means$Attention)
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)
# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion
dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))
# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"
# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
"Acq" = "Training",
"Bsln" = "Baseline",
"Ext" = "Test")
dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
"High_Rew" = "High reward probability",
"Low_Rew" = "Low reward probability")
dataPlot$Attention = recode(dataPlot$Attention,
"Att" = "Target",
"NotAtt" = "Distractor")
##### Plotting targets and distractors separately #####
# Attended
dataAmplitudePlot=subset(dataPlot,Attention=="Target")
posterior_means_plot=subset(posterior_means,Attention=="Target")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Target")
# tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Unattended
dataAmplitudePlot=subset(dataPlot,Attention=="Distractor")
posterior_means_plot=subset(posterior_means,Attention=="Distractor")
posterior_conditions_plot=subset(posterior_conditions,Attention=="Distractor")
# tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
##### Plotting high and low separately #####
# High
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="High reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
# Low
dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) +
# violin of the raw data
geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
labs(title="Low reward") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()
####### Plotting everything together #####
# Order for plotting
# dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor"))
# posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor"))
# posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor"))
#
# dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]
dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")
# tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)
ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(`Reward Probability`~Attention) +
scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
theme(
axis.line = element_line(size=2, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 24, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 22),
axis.text.x=element_text(colour="black", size = 22),
axis.text.y=element_text(colour="black", size = 22),
axis.title=element_text(size=24,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 24))
# dev.off()Table of means across conditions
# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward",
"Attended Baseline Low Reward",
"Attended Acquisition High Reward",
"Attended Acquisition Low Reward",
"Attended Extinction High Reward",
"Attended Extinction Low Reward",
"Unattended Baseline High Reward",
"Unattended Baseline Low Reward",
"Unattended Acquisition High Reward",
"Unattended Acquisition Low Reward",
"Unattended Extinction High Reward",
"Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Attention | Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|---|
| Attended | Baseline | High Reward | 1.14 [ 1.09 1.18 ] |
| Attended | Baseline | Low Reward | 1.13 [ 1.07 1.18 ] |
| Attended | Acquisition | High Reward | 1.17 [ 1.1 1.23 ] |
| Attended | Acquisition | Low Reward | 1.15 [ 1.08 1.22 ] |
| Attended | Extinction | High Reward | 1.13 [ 1.07 1.19 ] |
| Attended | Extinction | Low Reward | 0.91 [ 0.84 0.98 ] |
| Unattended | Baseline | High Reward | 0.89 [ 0.83 0.95 ] |
| Unattended | Baseline | Low Reward | 0.88 [ 0.82 0.94 ] |
| Unattended | Acquisition | High Reward | 0.92 [ 0.85 0.99 ] |
| Unattended | Acquisition | Low Reward | 0.9 [ 0.83 0.97 ] |
| Unattended | Extinction | High Reward | 0.88 [ 0.82 0.95 ] |
| Unattended | Extinction | Low Reward | 0.91 [ 0.84 0.98 ] |
| ### Inference | about the best | model |
Attended vs. unattended
Check the difference between attended and not attended in baseline high rewarded
Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in baseline low rewarded
Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition high rewarded
Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in acquisition low rewarded
Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction high rewarded
Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between attended and not attended in extinction low rewarded
Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Comparison between phases
Check the difference between baseline and acquisition in high reward attended
Diff_Bsln_Acq_High_Att = Baseline_High_Attended - Acquisition_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward attended
Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward not attended
Diff_Bsln_Acq_High_NotAtt = Baseline_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in low reward not attended
Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward attended
Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward attended
Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in high reward not attended
Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between acquisition and extinction in low reward not attended
Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Baseline difference
Check the difference between high and low reward in baseline attended
Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between high and low reward in baseline not attended
Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)
Check the difference between baseline and acquisition in high reward unattended vs. attended
Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")## [1] "Mean = 0 [ 0 0 ]"Behavior
Import data
# Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# # clear the environment
# rm(list=ls())
# # clear the console
# cat("\014")
# #load packages and install them if they're not installed
# if (!require("pacman")) install.packages("pacman")
# pacman::p_load(plyr,Rmisc,yarrr,BayesFactor,reshape2,brms, broom, tidyverse, brmstools, BEST, knitr, here,psych)
# # set seed
# set.seed(42)
# # set directory
# setwd(here())
# import data
data.raw = read.csv(file = here("data","Data_behavior_exp1_48pps.csv"),header=TRUE,na.strings="NaN")
# Prepare the dataset------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
### Adding and renaming variables
# rename EventType variable
names(data.raw)[names(data.raw) == "EventType"] = "MovedDots"
# add a variable with the name of the attended color instead of a numbers
data.raw$AttendedColor = ifelse(data.raw$AttendedColor==1,"red","blue")
# add a variable saying which color was linked with High_Rew (even numbers - blue was High_Rew)
data.raw$RewardedColor = ifelse(data.raw$ParticipantNo%%2==0,"blue","red")
# add a variable with the name of the moved color instead of a numbers
data.raw$MovedDots = ifelse(data.raw$MovedDots==1,"red","blue")
# split experimental phases into 6 isntead of 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
#data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,100,200,300,400,500,600),labels=c("Bsln1","Bsln2","Acq1","Acq2","Ext1","Ext2"))
# split experimental phases into 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,200,400,600),labels=c("Bsln","Acq","Ext")) # trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext
### Convert variables to be used in analyses into factors
data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )] =
lapply(data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )], factor)
### Create variables needed for the accuracy analyses
# count hits, false alarms, misses, correct rejections, and RT separately for each participant (their calculation is done in Matlab: see DataProcessing.m)
data.final = ddply(data.raw,.(ParticipantNo,ExpPhase,AttendedColor,RewardedColor,MovedDots),summarize,
numtrials=length(which(Response!=99)), # number of trials per condition (anything that is not 99 or any other number that we're not using)
Hits=length(which(Response==1)), # hits: attended color moved, correct response
FAs=length(which(Response==2)), # false alarms: attended color did not move, (wrong) response
Misses=length(which(Response==0)), # misses: attended color moved, no response
CRs=length(which(Response==3)), # correct rejections: attended color did not move, no response
mean.RT=mean(RT,na.rm=TRUE)) # mean RT per condition
################################################################## Calculate accuracy and RTs per condition ###############################################################################################################################################################################################################
# Prepare the data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
### Calculate Hits and False alarms
# Hits are calculated for each participant in each condition on trials when they are attending the color that moved.
# False alarms are calculated for each participant in each condition on trials when they are attending the color that didn't move (the unattended color moved, but they responded)
# Here we create the same number of hits & fas for each of the two conditions (moved attended or not)
data.final = ddply(data.final, .(ParticipantNo,ExpPhase,AttendedColor), transform,
Hits = Hits[MovedDots==AttendedColor],
FAs = FAs[MovedDots!=AttendedColor])
# Keep only trials on which the attended color moved (we can do behavioral analysis only on those)
data.final = subset(data.final,MovedDots==AttendedColor)
### Calculate d'
# use loglinear transformation: add 0.5 to Hits, FAs, Misses, and CRs (Hautus, 1995, Behavior Research Methods, Instruments, & Computers),
# which is preferred over the 1/2N rule (Macmillan & Kaplan, 1985, Psychological Bulletin) because it results in less biased estimates of d'.
data.final = ddply(data.final,.(ParticipantNo,ExpPhase,RewardedColor,AttendedColor,numtrials),summarise,
tot.Hits=Hits+.5, # hits
tot.FAs=FAs+.5, # false alarms
tot.Misses=(numtrials-Hits)+.5, # misses
tot.CRs=(numtrials-FAs)+.5, # correct rejections
Hit.Rate=tot.Hits/(tot.Hits+tot.Misses), # hit rate
FA.Rate=tot.FAs/(tot.FAs+tot.CRs), # false alarm rate
#dprime=qnorm(Hit.Rate)-qnorm(FA.Rate),
Hits.RTs=mean(mean.RT,na.rm=TRUE)) # mean RTs
# Calculate SDT indices with psycho
indices = psycho::dprime(data.final$tot.Hits, data.final$tot.FAs, data.final$tot.Misses, data.final$tot.CRs)
data.final = cbind(data.final, indices)
### Create a final dataframe for accuracy and RTs analyses
# add a new variable specifying whether the participant is attending the high or Low_Rewed color
data.final$Condition = ifelse(data.final$RewardedColor==data.final$AttendedColor,"High_Rew","Low_Rew")
# make this variable a factor for further analyses
data.final$Condition = factor(data.final$Condition)
# Exclude subjects without EEG
missing_eeg = c(5,10,13,27)
data.final = data.final[!data.final$ParticipantNo %in% missing_eeg,]D prime
Means - D prime - raw data
summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(dprime, na.rm = TRUE), digits = 2), " [", round(hdi(dprime)[[1]], digits = 2), " ", round(hdi(dprime)[[2]], digits = 2), "]")))
names(summary) = c("Reward phase", "Reward probability", "Dprime")
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
kable(summary, caption = "Dprime per condition")| Reward phase | Reward probability | Dprime |
|---|---|---|
| Baseline | High | 1.65 [ -0.04 2.68 ] |
| Baseline | Low | 1.45 [ 0.04 2.3 ] |
| Acquisition | High | 1.69 [ -0.29 2.73 ] |
| Acquisition | Low | 1.61 [ 0.46 2.68 ] |
| Extinction | High | 1.62 [ -0.2 2.73 ] |
| Extinction | Low | 1.61 [ 0.74 2.88 ] |
Plot
# # Plot FA rates------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] <- "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] <- "Reward probability"
# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
# Pirate plot
pirateplot(formula=dprime ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=data.plot, # data frame
main="Dprime", # main title
ylim=c(-2,4), # y-axis: limits
ylab="Dprime", # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.fun = median,
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.full.Acc.dprime = readRDS("model.full.Acc.dprime.rds")
compare.waic.Acc.dprime = readRDS("compare.Acc.dprime.waic")
compare.Acc.dprime.waic.weights = readRDS("compare.Acc.dprime.waic.weights")
bR2.null.Acc.dprime = readRDS("bR2.null.Acc.dprime")
bR2.expphase.Acc.dprime = readRDS("bR2.expphase.Acc.dprime")
bR2.full.Acc.dprime = readRDS("bR2.full.Acc.dprime")Model comparison with WAIC
print(compare.waic.Acc.dprime)## WAIC SE
## model.null.Acc.dprime 446.33 23.43
## model.expphase.Acc.dprime 451.47 22.74
## model.full.Acc.dprime 200.57 21.27
## model.null.Acc.dprime - model.expphase.Acc.dprime -5.14 3.35
## model.null.Acc.dprime - model.full.Acc.dprime 245.77 23.82
## model.expphase.Acc.dprime - model.full.Acc.dprime 250.91 23.53Model comparison with WAIC weights
print(compare.Acc.dprime.waic.weights)## model.null.Acc.dprime model.expphase.Acc.dprime
## 4.284133e-54 3.283575e-55
## model.full.Acc.dprime
## 1.000000e+00Bayesian R squared
Null model
print(bR2.null.Acc.dprime)## Estimate Est.Error Q2.5 Q97.5
## R2 0.2229344 0.04885937 0.125037 0.3154086Experiment phase model
print(bR2.expphase.Acc.dprime)## Estimate Est.Error Q2.5 Q97.5
## R2 0.2396632 0.04856887 0.1427286 0.3314708Full model
print(bR2.full.Acc.dprime)## Estimate Est.Error Q2.5 Q97.5
## R2 0.7620456 0.02220512 0.714308 0.8017077Checking the best model
Plotting the chains
# Plot chains
plot(model.full.Acc.dprime, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.full.Acc.dprime, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 1.74 0.09 1.592 1.88
## 2 ExpPhaseAcq 0.04 0.06 -0.063 0.15
## 3 ExpPhaseExt 0.01 0.06 -0.091 0.12
## 4 ConditionLow_Rew -0.24 0.14 -0.464 -0.02
## 5 ExpPhaseAcq:ConditionLow_Rew 0.14 0.09 -0.004 0.28
## 6 ExpPhaseExt:ConditionLow_Rew 0.17 0.09 0.024 0.32Plotting the best model
Sample from the posterior
# Analyzing the posterior and differences between conditions
post = posterior_samples(model.full.Acc.dprime, "^b")
################################################ Baseline ####
######### High reward
Baseline_High = post[["b_Intercept"]]
######### Low reward
Baseline_Low = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
################################################ Acquistion
######### High reward
Acquisition_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]]
######### Low reward
Acquisition_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq:ConditionLow_Rew"]]
################################################ Extinction
######### High reward
Extinction_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]]
######### Low reward
Extinction_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt:ConditionLow_Rew"]]# Data from the model
# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High, Baseline_Low, Acquisition_High, Acquisition_Low, Extinction_High, Extinction_Low))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")
names(posterior_conditions)[3] = "dprime"
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward",
"Low" = "Low reward")
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline High Reward",
"Baseline Low Reward",
"Acquisition High Reward",
"Acquisition Low Reward",
"Extinction High Reward",
"Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline_High)[[1]], digits = 2),
round(hdi(Baseline_Low)[[1]], digits = 2),
round(hdi(Acquisition_High)[[1]], digits = 2),
round(hdi(Acquisition_Low)[[1]], digits = 2),
round(hdi(Extinction_High)[[1]], digits = 2),
round(hdi(Extinction_Low)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline_High)[[2]], digits = 2),
round(hdi(Baseline_Low)[[2]], digits = 2),
round(hdi(Acquisition_High)[[2]], digits = 2),
round(hdi(Acquisition_Low)[[2]], digits = 2),
round(hdi(Extinction_High)[[2]], digits = 2),
round(hdi(Extinction_Low)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline_High,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction_High,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline_High,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction_High,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward",
"Low Reward" = "Low reward")
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`dprime` = c(1,1,1,1,1,1)
# Raw data
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"
# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
"Acq" = "Training",
"Bsln" = "Baseline",
"Ext" = "Test")
data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
"High_Rew" = "High reward",
"Low_Rew" = "Low reward")
# Violin + the model
# tiff("dprime.tiff", units="in", width=8, height=5, res=800)
ggplot(data=data.plot, aes(x=`Reward Phase`, y=`dprime`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`dprime`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_95_CI, ymax=`dprime`-`dprime`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_50_CI, ymax=`dprime`-`dprime`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(breaks=seq(-1,3,0.5)) +
labs(title="Sensitivity d`", y = "Sensitivity d`") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()Table of means across conditions
library(HDInterval)
# Make a table with conditions
posterior_means = as.data.frame(c("Baseline High Reward",
"Baseline Low Reward",
"Acquisition High Reward",
"Acquisition Low Reward",
"Extinction High Reward",
"Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline_High), digits = 2), " [", round(hdi(Baseline_High)[[1]], digits = 2), " ", round(hdi(Baseline_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low), digits = 2), " [", round(hdi(Baseline_Low)[[1]], digits = 2), " ", round(hdi(Baseline_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High), digits = 2), " [", round(hdi(Acquisition_High)[[1]], digits = 2), " ", round(hdi(Acquisition_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low), digits = 2), " [", round(hdi(Acquisition_Low)[[1]], digits = 2), " ", round(hdi(Acquisition_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High), digits = 2), " [", round(hdi(Extinction_High)[[1]], digits = 2), " ", round(hdi(Extinction_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low), digits = 2), " [", round(hdi(Extinction_Low)[[1]], digits = 2), " ", round(hdi(Extinction_Low)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|
| Baseline | High Reward | 1.74 [ 1.55 1.9 ] |
| Baseline | Low Reward | 1.5 [ 1.31 1.69 ] |
| Acquisition | High Reward | 1.78 [ 1.59 1.96 ] |
| Acquisition | Low Reward | 1.68 [ 1.49 1.87 ] |
| Extinction | High Reward | 1.75 [ 1.58 1.93 ] |
| Extinction | Low Reward | 1.68 [ 1.49 1.87 ] |
Inference about the best model
Check the difference between high reward in baseline vs. acquisition
Diff_Bsln_Acq_High = Acquisition_High - Baseline_High
plotPost(Diff_Bsln_Acq_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between low reward in baseline vs. acquisition
Diff_Bsln_Acq_Low = Acquisition_Low - Baseline_Low
plotPost(Diff_Bsln_Acq_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between high and low reward in baseline vs. acquisition
Diff_Bsln_Acq_High_vs_Low = Diff_Bsln_Acq_Low - Diff_Bsln_Acq_High
plotPost(Diff_Bsln_Acq_High_vs_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[2]], digits = 2), "]")## [1] "Mean = 0.14 [ -0.03 0.31 ]"Check the difference between high reward in acquisition vs. extinction
Diff_Acq_Ext_High = Extinction_High - Acquisition_High
plotPost(Diff_Acq_Ext_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between low reward in acquisition vs. extinction
Diff_Acq_Ext_Low = Extinction_Low - Acquisition_Low
plotPost(Diff_Acq_Ext_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Reaction times
Means - raw data
summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Hits.RTs, na.rm = TRUE), digits = 2), " [", round(hdi(Hits.RTs)[[1]], digits = 2), " ", round(hdi(Hits.RTs)[[2]], digits = 2), "]")))
names(summary) = c("Reward phase", "Reward probability", "Reaction times")
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
kable(summary, caption = "Reaction times per condition")| Reward phase | Reward probability | Reaction times |
|---|---|---|
| Baseline | High | 547.07 [ 485.64 619.34 ] |
| Baseline | Low | 553.88 [ 480.45 631.36 ] |
| Acquisition | High | 525.11 [ 467.12 599.49 ] |
| Acquisition | Low | 538.55 [ 465.32 584.63 ] |
| Extinction | High | 528.66 [ 457.08 599.83 ] |
| Extinction | Low | 539.76 [ 455.8 623.21 ] |
Plot
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward probability"
# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
# Pirate plot
pirateplot(formula = Hits.RTs ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=data.plot, # data frame
main="Reaction times", # main title
ylim=c(400,700), # y-axis: limits
ylab="Reaction time", # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.fun = median,
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.full.RT = readRDS("model.full.RT.rds")
compare.waic.RT = readRDS("compare.RT.waic")
compare.RT.waic.weights = readRDS("compare.RT.waic.weights")
bR2.null.RT = readRDS("bR2.null.RT")
bR2.expphase.RT = readRDS("bR2.expphase.RT")
bR2.full.RT = readRDS("bR2.full.RT")Model comparison with WAIC
print(compare.waic.RT)## WAIC SE
## model.null.RT 2541.30 26.80
## model.expphase.RT 2511.86 27.49
## model.full.RT 2324.94 23.46
## model.null.RT - model.expphase.RT 29.44 11.92
## model.null.RT - model.full.RT 216.36 18.93
## model.expphase.RT - model.full.RT 186.92 17.80Bayesian R squared
Null model
print(bR2.null.RT)## Estimate Est.Error Q2.5 Q97.5
## R2 0.584279 0.03049892 0.5184536 0.6374979Experiment phase model
print(bR2.expphase.RT)## Estimate Est.Error Q2.5 Q97.5
## R2 0.6484994 0.02806138 0.5888214 0.6997926Full model
print(bR2.full.RT)## Estimate Est.Error Q2.5 Q97.5
## R2 0.8666485 0.01580508 0.8315594 0.8934709Checking the best model
Plotting the chains
# Plot chains
plot(model.full.RT, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.full.RT, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 551 6 540.9 561
## 2 ExpPhaseAcq -22 4 -29.0 -15
## 3 ExpPhaseExt -18 5 -26.4 -10
## 4 ConditionLow_Rew 6 6 -4.0 15
## 5 ExpPhaseAcq:ConditionLow_Rew 7 5 -0.3 15
## 6 ExpPhaseExt:ConditionLow_Rew 2 5 -5.6 10Plotting the best model
Sample from the posterior
# Analyzing the posterior and differences between conditions
post = posterior_samples(model.full.RT, "^b")
################################################ Baseline ####
######### High reward
Baseline_High = post[["b_Intercept"]]
######### Low reward
Baseline_Low = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
################################################ Acquistion
######### High reward
Acquisition_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]]
######### Low reward
Acquisition_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq:ConditionLow_Rew"]]
################################################ Extinction
######### High reward
Extinction_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]]
######### Low reward
Extinction_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt:ConditionLow_Rew"]]# Data from the model
# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High, Baseline_Low, Acquisition_High, Acquisition_Low, Extinction_High, Extinction_Low))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")
names(posterior_conditions)[3] = "Reaction time"
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward",
"Low" = "Low reward")
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline High Reward",
"Baseline Low Reward",
"Acquisition High Reward",
"Acquisition Low Reward",
"Extinction High Reward",
"Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline_High)[[1]], digits = 2),
round(hdi(Baseline_Low)[[1]], digits = 2),
round(hdi(Acquisition_High)[[1]], digits = 2),
round(hdi(Acquisition_Low)[[1]], digits = 2),
round(hdi(Extinction_High)[[1]], digits = 2),
round(hdi(Extinction_Low)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline_High)[[2]], digits = 2),
round(hdi(Baseline_Low)[[2]], digits = 2),
round(hdi(Acquisition_High)[[2]], digits = 2),
round(hdi(Acquisition_Low)[[2]], digits = 2),
round(hdi(Extinction_High)[[2]], digits = 2),
round(hdi(Extinction_Low)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline_High,0.5)[[1]], digits = 2),
round(hdi(Baseline_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction_High,0.5)[[1]], digits = 2),
round(hdi(Extinction_Low,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline_High,0.5)[[2]], digits = 2),
round(hdi(Baseline_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction_High,0.5)[[2]], digits = 2),
round(hdi(Extinction_Low,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward",
"Low Reward" = "Low reward")
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition" = "Training",
"Baseline" = "Baseline",
"Extinction" = "Test")
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`Reaction time` = c(555,555,555,555,555,555)
# Raw data
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"
# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
"Acq" = "Training",
"Bsln" = "Baseline",
"Ext" = "Test")
data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
"High_Rew" = "High reward",
"Low_Rew" = "Low reward")
data.plot$`Reaction time` = data.plot$Hits.RTs
# Violin + the model
# tiff("RTs.tiff", units="in", width=8, height=5, res=800)
ggplot(data=data.plot, aes(x=`Reward Phase`, y=`Reaction time`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Reaction time`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_95_CI, ymax=`Reaction time`-`Reaction time`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_50_CI, ymax=`Reaction time`-`Reaction time`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(breaks=seq(400,750,50)) +
labs(title="Reaction time", y = "Reaction time") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()Table of means across conditions
library(HDInterval)
# Make a table with conditions
posterior_means = as.data.frame(c("Baseline High Reward",
"Baseline Low Reward",
"Acquisition High Reward",
"Acquisition Low Reward",
"Extinction High Reward",
"Extinction Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline_High), digits = 2), " [", round(hdi(Baseline_High)[[1]], digits = 2), " ", round(hdi(Baseline_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline_Low), digits = 2), " [", round(hdi(Baseline_Low)[[1]], digits = 2), " ", round(hdi(Baseline_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_High), digits = 2), " [", round(hdi(Acquisition_High)[[1]], digits = 2), " ", round(hdi(Acquisition_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition_Low), digits = 2), " [", round(hdi(Acquisition_Low)[[1]], digits = 2), " ", round(hdi(Acquisition_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_High), digits = 2), " [", round(hdi(Extinction_High)[[1]], digits = 2), " ", round(hdi(Extinction_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction_Low), digits = 2), " [", round(hdi(Extinction_Low)[[1]], digits = 2), " ", round(hdi(Extinction_Low)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|
| Baseline | High Reward | 551.15 [ 538.76 563.4 ] |
| Baseline | Low Reward | 556.82 [ 543.65 570.1 ] |
| Acquisition | High Reward | 528.87 [ 516.89 541.25 ] |
| Acquisition | Low Reward | 542.01 [ 530.03 554.53 ] |
| Extinction | High Reward | 533.06 [ 518.87 547.36 ] |
| Extinction | Low Reward | 540.89 [ 526.07 555.71 ] |
Inference about the best model
Check the difference between high reward in baseline vs. acquisition
Diff_Bsln_Acq_High = Acquisition_High - Baseline_High
plotPost(Diff_Bsln_Acq_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between low reward in baseline vs. acquisition
Diff_Bsln_Acq_Low = Acquisition_Low - Baseline_Low
plotPost(Diff_Bsln_Acq_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between high and low reward in baseline vs. acquisition
Diff_Bsln_Acq_High_vs_Low = Diff_Bsln_Acq_High - Diff_Bsln_Acq_Low
plotPost(Diff_Bsln_Acq_High_vs_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[2]], digits = 2), "]")## [1] "Mean = -7.46 [ -16.49 1.99 ]"Check the difference between high reward in acquisition vs. extinction
Diff_Acq_Ext_High = Extinction_High - Acquisition_High
plotPost(Diff_Acq_Ext_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between low reward in acquisition vs. extinction
Diff_Acq_Ext_Low = Extinction_Low - Acquisition_Low
plotPost(Diff_Acq_Ext_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Behavior - splitting the phases
Import data
# Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# # clear the environment
# rm(list=ls())
# # clear the console
# cat("\014")
# #load packages and install them if they're not installed
# if (!require("pacman")) install.packages("pacman")
# pacman::p_load(plyr,Rmisc,yarrr,BayesFactor,reshape2,brms, broom, tidyverse, brmstools, BEST, knitr, here,psych)
# # set seed
# set.seed(42)
# # set directory
# setwd(here())
# import data
data.raw = read.csv(file = here("data","Data_behavior_exp1_48pps.csv"),header=TRUE,na.strings="NaN")
# Prepare the dataset------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
### Adding and renaming variables
# rename EventType variable
names(data.raw)[names(data.raw) == "EventType"] = "MovedDots"
# add a variable with the name of the attended color instead of a numbers
data.raw$AttendedColor = ifelse(data.raw$AttendedColor==1,"red","blue")
# add a variable saying which color was linked with High_Rew (even numbers - blue was High_Rew)
data.raw$RewardedColor = ifelse(data.raw$ParticipantNo%%2==0,"blue","red")
# add a variable with the name of the moved color instead of a numbers
data.raw$MovedDots = ifelse(data.raw$MovedDots==1,"red","blue")
# split experimental phases into 6 isntead of 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,100,200,300,400,500,600),labels=c("Bsln1","Bsln2","Acq1","Acq2","Ext1","Ext2"))
# split experimental phases into 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
# data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,200,400,600),labels=c("Bsln","Acq","Ext")) # trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext
### Convert variables to be used in analyses into factors
data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )] =
lapply(data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )], factor)
### Create variables needed for the accuracy analyses
# count hits, false alarms, misses, correct rejections, and RT separately for each participant (their calculation is done in Matlab: see DataProcessing.m)
data.final = ddply(data.raw,.(ParticipantNo,ExpPhase,AttendedColor,RewardedColor,MovedDots),summarize,
numtrials=length(which(Response!=99)), # number of trials per condition (anything that is not 99 or any other number that we're not using)
Hits=length(which(Response==1)), # hits: attended color moved, correct response
FAs=length(which(Response==2)), # false alarms: attended color did not move, (wrong) response
Misses=length(which(Response==0)), # misses: attended color moved, no response
CRs=length(which(Response==3)), # correct rejections: attended color did not move, no response
mean.RT=mean(RT,na.rm=TRUE)) # mean RT per condition
################################################################## Calculate accuracy and RTs per condition ###############################################################################################################################################################################################################
# Prepare the data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
### Calculate Hits and False alarms
# Hits are calculated for each participant in each condition on trials when they are attending the color that moved.
# False alarms are calculated for each participant in each condition on trials when they are attending the color that didn't move (the unattended color moved, but they responded)
# Here we create the same number of hits & fas for each of the two conditions (moved attended or not)
data.final = ddply(data.final, .(ParticipantNo,ExpPhase,AttendedColor), transform,
Hits = Hits[MovedDots==AttendedColor],
FAs = FAs[MovedDots!=AttendedColor])
# Keep only trials on which the attended color moved (we can do behavioral analysis only on those)
data.final = subset(data.final,MovedDots==AttendedColor)
### Calculate d'
# use loglinear transformation: add 0.5 to Hits, FAs, Misses, and CRs (Hautus, 1995, Behavior Research Methods, Instruments, & Computers),
# which is preferred over the 1/2N rule (Macmillan & Kaplan, 1985, Psychological Bulletin) because it results in less biased estimates of d'.
data.final = ddply(data.final,.(ParticipantNo,ExpPhase,RewardedColor,AttendedColor,numtrials),summarise,
tot.Hits=Hits+.5, # hits
tot.FAs=FAs+.5, # false alarms
tot.Misses=(numtrials-Hits)+.5, # misses
tot.CRs=(numtrials-FAs)+.5, # correct rejections
Hit.Rate=tot.Hits/(tot.Hits+tot.Misses), # hit rate
FA.Rate=tot.FAs/(tot.FAs+tot.CRs), # false alarm rate
#dprime=qnorm(Hit.Rate)-qnorm(FA.Rate),
Hits.RTs=mean(mean.RT,na.rm=TRUE)) # mean RTs
# Calculate SDT indices with psycho
indices = psycho::dprime(data.final$tot.Hits, data.final$tot.FAs, data.final$tot.Misses, data.final$tot.CRs)
data.final = cbind(data.final, indices)
### Create a final dataframe for accuracy and RTs analyses
# add a new variable specifying whether the participant is attending the high or Low_Rewed color
data.final$Condition = ifelse(data.final$RewardedColor==data.final$AttendedColor,"High_Rew","Low_Rew")
# make this variable a factor for further analyses
data.final$Condition = factor(data.final$Condition)
# Exclude subjects without EEG
missing_eeg = c(5,10,13,27)
data.final = data.final[!data.final$ParticipantNo %in% missing_eeg,]D prime
Means - D prime - raw data
summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(dprime, na.rm = TRUE), digits = 2), " [", round(hdi(dprime)[[1]], digits = 2), " ", round(hdi(dprime)[[2]], digits = 2), "]")))
names(summary) = c("Reward phase", "Reward probability", "Dprime")
# rename conditions
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq1" = "Acquisition1",
"Acq2" = "Acquisition2",
"Bsln1" = "Baseline1",
"Bsln2" = "Baseline2",
"Ext1" = "Extinction1",
"Ext2" = "Extinction2")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
kable(summary, caption = "Dprime per condition")| Reward phase | Reward probability | Dprime |
|---|---|---|
| Baseline1 | High | 1.49 [ -0.36 2.62 ] |
| Baseline1 | Low | 1.3 [ 0.09 2.35 ] |
| Baseline2 | High | 1.6 [ -0.27 2.56 ] |
| Baseline2 | Low | 1.45 [ 0.08 2.33 ] |
| Acquisition1 | High | 1.56 [ -0.08 2.65 ] |
| Acquisition1 | Low | 1.58 [ 0.47 2.45 ] |
| Acquisition2 | High | 1.59 [ 0.08 2.56 ] |
| Acquisition2 | Low | 1.47 [ 0 2.62 ] |
| Extinction1 | High | 1.51 [ -0.07 2.57 ] |
| Extinction1 | Low | 1.49 [ 0.36 2.5 ] |
| Extinction2 | High | 1.5 [ -0.38 2.49 ] |
| Extinction2 | Low | 1.53 [ 0.65 2.55 ] |
Plot
# # Plot FA rates------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] <- "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] <- "Reward probability"
# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
"Acq1" = "Acquisition1",
"Acq2" = "Acquisition2",
"Bsln1" = "Baseline1",
"Bsln2" = "Baseline2",
"Ext1" = "Extinction1",
"Ext2" = "Extinction2")
data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
# Pirate plot
pirateplot(formula=dprime ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=data.plot, # data frame
main="Dprime", # main title
ylim=c(-2,4), # y-axis: limits
ylab="Dprime", # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.fun = median,
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.full.Acc.dprime.split = readRDS("model.full.Acc.dprime.split.rds")
compare.waic.Acc.dprime.split = readRDS("compare.Acc.dprime.split.waic")
compare.Acc.dprime.split.waic.weights = readRDS("compare.Acc.dprime.split.waic.weights")
bR2.null.Acc.dprime.split = readRDS("bR2.null.Acc.dprime.split")
bR2.expphase.Acc.dprime.split = readRDS("bR2.expphase.Acc.dprime.split")
bR2.full.Acc.dprime.split = readRDS("bR2.full.Acc.dprime.split")Model comparison with WAIC
print(compare.waic.Acc.dprime.split)## WAIC SE
## model.null.Acc.dprime 1027.66 35.66
## model.expphase.Acc.dprime 1044.58 35.48
## model.full.Acc.dprime 568.84 33.83
## model.null.Acc.dprime - model.expphase.Acc.dprime -16.92 4.73
## model.null.Acc.dprime - model.full.Acc.dprime 458.82 39.60
## model.expphase.Acc.dprime - model.full.Acc.dprime 475.74 39.58Model comparison with WAIC weights
print(compare.Acc.dprime.split.waic.weights)## model.null.Acc.dprime model.expphase.Acc.dprime
## 2.337776e-100 4.953552e-104
## model.full.Acc.dprime
## 1.000000e+00Bayesian R squared
Null model
print(bR2.null.Acc.dprime.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.2793751 0.03088977 0.2176021 0.3388459Experiment phase model
print(bR2.expphase.Acc.dprime.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.2943567 0.03005744 0.2333531 0.351622Full model
print(bR2.full.Acc.dprime.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.745154 0.01374807 0.7164296 0.7701522Checking the best model
Plotting the chains
# Plot chains
plot(model.full.Acc.dprime.split, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.full.Acc.dprime.split, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 1.49 0.12 1.29 1.68
## 2 ExpPhaseBsln2 0.12 0.08 -0.02 0.25
## 3 ExpPhaseAcq1 0.07 0.08 -0.06 0.21
## 4 ExpPhaseAcq2 0.10 0.08 -0.04 0.23
## 5 ExpPhaseExt1 0.02 0.08 -0.11 0.15
## 6 ExpPhaseExt2 0.01 0.08 -0.12 0.14
## 7 ConditionLow_Rew -0.19 0.16 -0.45 0.08
## 8 ExpPhaseBsln2:ConditionLow_Rew 0.03 0.11 -0.15 0.22
## 9 ExpPhaseAcq1:ConditionLow_Rew 0.21 0.11 0.03 0.39
## 10 ExpPhaseAcq2:ConditionLow_Rew 0.07 0.11 -0.11 0.26
## 11 ExpPhaseExt1:ConditionLow_Rew 0.17 0.11 -0.01 0.36
## 12 ExpPhaseExt2:ConditionLow_Rew 0.23 0.11 0.04 0.41Plotting the best model
Sample from the posterior
# Analyzing the posterior and differences between conditions
post = posterior_samples(model.full.Acc.dprime.split, "^b")
################################################ Baseline 1 ####
######### High reward
Baseline1_High = post[["b_Intercept"]]
######### Low reward
Baseline1_Low = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
################################################ Baseline 2 ####
######### High reward
Baseline2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseBsln2"]]
######### Low reward
Baseline2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseBsln2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseBsln2:ConditionLow_Rew"]]
################################################ Acquistion 1
######### High reward
Acquisition1_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq1"]]
######### Low reward
Acquisition1_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq1"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq1:ConditionLow_Rew"]]
################################################ Acquistion 2
######### High reward
Acquisition2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq2"]]
######### Low reward
Acquisition2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq2:ConditionLow_Rew"]]
################################################ Extinction 1
######### High reward
Extinction1_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt1"]]
######### Low reward
Extinction1_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt1"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt1:ConditionLow_Rew"]]
################################################ Extinction 2
######### High reward
Extinction2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt2"]]
######### Low reward
Extinction2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt2:ConditionLow_Rew"]]# Data from the model
# make a data frame
posterior_conditions = melt(data.frame(Baseline1_High, Baseline2_High, Baseline1_Low, Baseline2_Low, Acquisition1_High, Acquisition1_Low, Acquisition2_High, Acquisition2_Low, Extinction1_High, Extinction1_Low, Extinction2_High, Extinction2_Low))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")
names(posterior_conditions)[3] = "dprime"
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward",
"Low" = "Low reward")
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition1" = "Training1",
"Acquisition2" = "Training2",
"Baseline1" = "Baseline1",
"Baseline1" = "Baseline2",
"Extinction1" = "Test1",
"Extinction2" = "Test2")
# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline1 High Reward",
"Baseline2 High Reward",
"Baseline1 Low Reward",
"Baseline2 Low Reward",
"Acquisition1 High Reward",
"Acquisition2 High Reward",
"Acquisition1 Low Reward",
"Acquisition2 Low Reward",
"Extinction1 High Reward",
"Extinction2 High Reward",
"Extinction1 Low Reward",
"Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline1_High)[[1]], digits = 2),
round(hdi(Baseline2_High)[[1]], digits = 2),
round(hdi(Baseline1_Low)[[1]], digits = 2),
round(hdi(Baseline2_Low)[[1]], digits = 2),
round(hdi(Acquisition1_High)[[1]], digits = 2),
round(hdi(Acquisition2_High)[[1]], digits = 2),
round(hdi(Acquisition1_Low)[[1]], digits = 2),
round(hdi(Acquisition2_Low)[[1]], digits = 2),
round(hdi(Extinction1_High)[[1]], digits = 2),
round(hdi(Extinction2_High)[[1]], digits = 2),
round(hdi(Extinction1_Low)[[1]], digits = 2),
round(hdi(Extinction2_Low)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline1_High)[[2]], digits = 2),
round(hdi(Baseline2_High)[[2]], digits = 2),
round(hdi(Baseline1_Low)[[2]], digits = 2),
round(hdi(Baseline2_Low)[[2]], digits = 2),
round(hdi(Acquisition1_High)[[2]], digits = 2),
round(hdi(Acquisition2_High)[[2]], digits = 2),
round(hdi(Acquisition1_Low)[[2]], digits = 2),
round(hdi(Acquisition2_Low)[[2]], digits = 2),
round(hdi(Extinction1_High)[[2]], digits = 2),
round(hdi(Extinction2_High)[[2]], digits = 2),
round(hdi(Extinction1_Low)[[2]], digits = 2),
round(hdi(Extinction2_Low)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline1_High,0.5)[[1]], digits = 2),
round(hdi(Baseline2_High,0.5)[[1]], digits = 2),
round(hdi(Baseline1_Low,0.5)[[1]], digits = 2),
round(hdi(Baseline2_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition1_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition2_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition1_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition2_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction1_High,0.5)[[1]], digits = 2),
round(hdi(Extinction2_High,0.5)[[1]], digits = 2),
round(hdi(Extinction1_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction2_Low,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline1_High,0.5)[[2]], digits = 2),
round(hdi(Baseline2_High,0.5)[[2]], digits = 2),
round(hdi(Baseline1_Low,0.5)[[2]], digits = 2),
round(hdi(Baseline2_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition1_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition2_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition1_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition2_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction1_High,0.5)[[2]], digits = 2),
round(hdi(Extinction2_High,0.5)[[2]], digits = 2),
round(hdi(Extinction1_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction2_Low,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward",
"Low Reward" = "Low reward")
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition1" = "Training1",
"Acquisition2" = "Training2",
"Baseline1" = "Baseline1",
"Baseline1" = "Baseline2",
"Extinction1" = "Test1",
"Extinction2" = "Test2")
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`dprime` = c(1,1,1,1,1,1)
# Raw data
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"
# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
"Acq1" = "Training1",
"Acq2" = "Training2",
"Bsln1" = "Baseline1",
"Bsln2" = "Baseline2",
"Ext1" = "Test1",
"Ext2" = "Test2")
data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
"High_Rew" = "High reward",
"Low_Rew" = "Low reward")
# Violin + the model
# tiff("dprime_split.tiff", units="in", width=16, height=5, res=800)
ggplot(data=data.plot, aes(x=`Reward Phase`, y=`dprime`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`dprime`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_95_CI, ymax=`dprime`-`dprime`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_50_CI, ymax=`dprime`-`dprime`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline1","Baseline2", "Training1","Training2", "Test1","Test2")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(breaks=seq(-1,3,0.5)) +
labs(title="Sensitivity d`", y = "Sensitivity d`") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()Table of means across conditions
# Make a table with conditions
posterior_means = as.data.frame(c("Baseline1 High Reward",
"Baseline1 Low Reward",
"Baseline2 High Reward",
"Baseline2 Low Reward",
"Acquisition1 High Reward",
"Acquisition1 Low Reward",
"Acquisition2 High Reward",
"Acquisition2 Low Reward",
"Extinction1 High Reward",
"Extinction1 Low Reward",
"Extinction2 High Reward",
"Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline1_High), digits = 2), " [", round(hdi(Baseline1_High)[[1]], digits = 2), " ", round(hdi(Baseline1_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline1_Low), digits = 2), " [", round(hdi(Baseline1_Low)[[1]], digits = 2), " ", round(hdi(Baseline1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Baseline2_High), digits = 2), " [", round(hdi(Baseline2_High)[[1]], digits = 2), " ", round(hdi(Baseline2_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline2_Low), digits = 2), " [", round(hdi(Baseline2_Low)[[1]], digits = 2), " ", round(hdi(Baseline2_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition1_High), digits = 2), " [", round(hdi(Acquisition1_High)[[1]], digits = 2), " ", round(hdi(Acquisition1_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition1_Low), digits = 2), " [", round(hdi(Acquisition1_Low)[[1]], digits = 2), " ", round(hdi(Acquisition1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition2_High), digits = 2), " [", round(hdi(Acquisition2_High)[[1]], digits = 2), " ", round(hdi(Acquisition2_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition2_Low), digits = 2), " [", round(hdi(Acquisition2_Low)[[1]], digits = 2), " ", round(hdi(Acquisition2_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction1_High), digits = 2), " [", round(hdi(Extinction1_High)[[1]], digits = 2), " ", round(hdi(Extinction1_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction1_Low), digits = 2), " [", round(hdi(Extinction1_Low)[[1]], digits = 2), " ", round(hdi(Extinction1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction2_High), digits = 2), " [", round(hdi(Extinction2_High)[[1]], digits = 2), " ", round(hdi(Extinction2_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction2_Low), digits = 2), " [", round(hdi(Extinction2_Low)[[1]], digits = 2), " ", round(hdi(Extinction2_Low)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|
| Baseline1 | High Reward | 1.49 [ 1.26 1.71 ] |
| Baseline1 | Low Reward | 1.3 [ 1.08 1.51 ] |
| Baseline2 | High Reward | 1.6 [ 1.38 1.83 ] |
| Baseline2 | Low Reward | 1.45 [ 1.23 1.66 ] |
| Acquisition1 | High Reward | 1.56 [ 1.32 1.79 ] |
| Acquisition1 | Low Reward | 1.58 [ 1.38 1.8 ] |
| Acquisition2 | High Reward | 1.59 [ 1.34 1.81 ] |
| Acquisition2 | Low Reward | 1.47 [ 1.24 1.68 ] |
| Extinction1 | High Reward | 1.51 [ 1.28 1.74 ] |
| Extinction1 | Low Reward | 1.49 [ 1.28 1.71 ] |
| Extinction2 | High Reward | 1.5 [ 1.25 1.72 ] |
| Extinction2 | Low Reward | 1.54 [ 1.32 1.75 ] |
Inference about the best model
Check the difference between baseline1 vs. baseline2 in high reward
Diff_Bsln1_Bsln2_High = Baseline2_High - Baseline1_High
plotPost(Diff_Bsln1_Bsln2_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline1 vs. baseline2 in low reward
Diff_Bsln1_Bsln2_Low = Baseline2_Low - Baseline1_Low
plotPost(Diff_Bsln1_Bsln2_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline2 vs. acquisition1 in high reward
Diff_Bsln2_Acq1_High = Baseline2_High - Acquisition1_High
plotPost(Diff_Bsln2_Acq1_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline2 vs. acquisition1 in low reward
Diff_Bsln2_Acq1_Low = Acquisition1_Low - Baseline2_Low
plotPost(Diff_Bsln2_Acq1_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Reaction times
Means - raw data
summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Hits.RTs, na.rm = TRUE), digits = 2), " [", round(hdi(Hits.RTs)[[1]], digits = 2), " ", round(hdi(Hits.RTs)[[2]], digits = 2), "]")))
names(summary) = c("Reward phase", "Reward probability", "Reaction times")
# rename conditions
summary$`Reward phase` = recode(summary$`Reward phase`,
"Acq1" = "Acquisition1",
"Acq2" = "Acquisition2",
"Bsln1" = "Baseline1",
"Bsln2" = "Baseline2",
"Ext1" = "Extinction1",
"Ext2" = "Extinction2")
summary$`Reward probability` = recode(summary$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
summary = as.data.frame(summary)
kable(summary, caption = "Reaction times per condition")| Reward phase | Reward probability | Reaction times |
|---|---|---|
| Baseline1 | High | 548.73 [ 479.43 613.76 ] |
| Baseline1 | Low | 551.97 [ 458.26 621.44 ] |
| Baseline2 | High | 545.6 [ 454.56 620.36 ] |
| Baseline2 | Low | 556.21 [ 486.84 650.73 ] |
| Acquisition1 | High | 521.47 [ 437.9 587.57 ] |
| Acquisition1 | Low | 542.8 [ 463.65 593.47 ] |
| Acquisition2 | High | 528.96 [ 462 598.58 ] |
| Acquisition2 | Low | 534.72 [ 479.38 618.25 ] |
| Extinction1 | High | 528.3 [ 457.88 596.17 ] |
| Extinction1 | Low | 536.86 [ 444.89 621 ] |
| Extinction2 | High | 528.92 [ 448.24 606 ] |
| Extinction2 | Low | 542.76 [ 450.11 617.44 ] |
Plot
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward probability"
# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
"Acq" = "Acquisition",
"Bsln" = "Baseline",
"Ext" = "Extinction")
data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
"High_Rew" = "High",
"Low_Rew" = "Low")
# Pirate plot
pirateplot(formula = Hits.RTs ~ `Reward phase` + `Reward probability`, # dependent~independent variables
data=data.plot, # data frame
main="Reaction times", # main title
ylim=c(400,700), # y-axis: limits
ylab="Reaction time", # y-axis: label
theme=0, # preset theme (0: use your own)
point.col="black", # points: color
point.o=.3, # points: opacity (0-1)
avg.line.fun = median,
avg.line.col="black", # average line: color
avg.line.lwd=2, # average line: line width
avg.line.o=1, # average line: opacity (0-1)
bean.b.col="black", # bean border, color
bean.lwd=0.6, # bean border, line width
bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
bean.b.o=0.3, # bean border, opacity (0-1)
bean.f.col="gray", # bean filling, color
bean.f.o=.1, # bean filling, opacity (0-1)
cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
gl.col="gray", # gridlines: color
gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
cex.lab=1, # axis labels: size
cex.axis=1, # axis numbers: size
cex.names = 1,
bty="l", # plot box type
back.col="white") # background, color
Statistics
# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))
model.full.RT.split = readRDS("model.full.RT.split.rds")
compare.waic.RT.split = readRDS("compare.RT.split.waic")
compare.RT.split.waic.weights = readRDS("compare.RT.split.waic.weights")
bR2.null.RT.split = readRDS("bR2.null.RT.split")
bR2.expphase.RT.split = readRDS("bR2.expphase.RT.split")
bR2.full.RT.split = readRDS("bR2.full.RT.split")Model comparison with WAIC
print(compare.waic.RT.split)## WAIC SE
## model.null.RT 5222.80 50.62
## model.expphase.RT 5191.25 53.84
## model.full.RT 4985.22 52.53
## model.null.RT - model.expphase.RT 31.55 16.03
## model.null.RT - model.full.RT 237.58 27.39
## model.expphase.RT - model.full.RT 206.03 25.45Bayesian R squared
Null model
print(bR2.null.RT.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.4448848 0.02600351 0.3915422 0.4929036Experiment phase model
print(bR2.expphase.RT.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.5320485 0.02751194 0.4752605 0.5816757Full model
print(bR2.full.RT.split)## Estimate Est.Error Q2.5 Q97.5
## R2 0.7264248 0.0185819 0.6871889 0.7599188Checking the best model
Plotting the chains
# Plot chains
plot(model.full.RT.split, pars = "^b_", ask = FALSE, N=6)
Summary of the best model
# Summary of the best model
print(tidy(model.full.RT.split, par_type = "non-varying"), digits = 1)## term estimate std.error lower upper
## 1 Intercept 549 6 538.2 559
## 2 ExpPhaseBsln2 -3 5 -11.7 6
## 3 ExpPhaseAcq1 -27 5 -35.5 -19
## 4 ExpPhaseAcq2 -20 5 -28.5 -11
## 5 ExpPhaseExt1 -20 6 -30.1 -11
## 6 ExpPhaseExt2 -20 6 -29.5 -10
## 7 ConditionLow_Rew 3 7 -7.9 15
## 8 ExpPhaseBsln2:ConditionLow_Rew 7 7 -4.1 19
## 9 ExpPhaseAcq1:ConditionLow_Rew 18 7 6.4 29
## 10 ExpPhaseAcq2:ConditionLow_Rew 2 7 -9.1 14
## 11 ExpPhaseExt1:ConditionLow_Rew 5 7 -6.3 17
## 12 ExpPhaseExt2:ConditionLow_Rew 11 7 -0.7 22Plotting the best model
Sample from the posterior
# Analyzing the posterior and differences between conditions
post = posterior_samples(model.full.RT.split, "^b")
################################################ Baseline 1 ####
######### High reward
Baseline1_High = post[["b_Intercept"]]
######### Low reward
Baseline1_Low = post[["b_Intercept"]] +
post[["b_ConditionLow_Rew"]]
################################################ Baseline 2 ####
######### High reward
Baseline2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseBsln2"]]
######### Low reward
Baseline2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseBsln2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseBsln2:ConditionLow_Rew"]]
################################################ Acquistion 1
######### High reward
Acquisition1_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq1"]]
######### Low reward
Acquisition1_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq1"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq1:ConditionLow_Rew"]]
################################################ Acquistion 2
######### High reward
Acquisition2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq2"]]
######### Low reward
Acquisition2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseAcq2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseAcq2:ConditionLow_Rew"]]
################################################ Extinction 1
######### High reward
Extinction1_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt1"]]
######### Low reward
Extinction1_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt1"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt1:ConditionLow_Rew"]]
################################################ Extinction 2
######### High reward
Extinction2_High = post[["b_Intercept"]] +
post[["b_ExpPhaseExt2"]]
######### Low reward
Extinction2_Low = post[["b_Intercept"]] +
post[["b_ExpPhaseExt2"]] +
post[["b_ConditionLow_Rew"]] +
post[["b_ExpPhaseExt2:ConditionLow_Rew"]]# Data from the model
# make a data frame
posterior_conditions = melt(data.frame(Baseline1_High, Baseline2_High, Baseline1_Low, Baseline2_Low, Acquisition1_High, Acquisition1_Low, Acquisition2_High, Acquisition2_Low, Extinction1_High, Extinction1_Low, Extinction2_High, Extinction2_Low))
posterior_conditions = posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")
names(posterior_conditions)[3] = "Reaction time"
posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
"High" = "High reward",
"Low" = "Low reward")
posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
"Acquisition1" = "Training1",
"Acquisition2" = "Training2",
"Baseline1" = "Baseline1",
"Baseline1" = "Baseline2",
"Extinction1" = "Test1",
"Extinction2" = "Test2")
# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline1 High Reward",
"Baseline2 High Reward",
"Baseline1 Low Reward",
"Baseline2 Low Reward",
"Acquisition1 High Reward",
"Acquisition2 High Reward",
"Acquisition1 Low Reward",
"Acquisition2 Low Reward",
"Extinction1 High Reward",
"Extinction2 High Reward",
"Extinction1 Low Reward",
"Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$low_95_CI = c(round(hdi(Baseline1_High)[[1]], digits = 2),
round(hdi(Baseline2_High)[[1]], digits = 2),
round(hdi(Baseline1_Low)[[1]], digits = 2),
round(hdi(Baseline2_Low)[[1]], digits = 2),
round(hdi(Acquisition1_High)[[1]], digits = 2),
round(hdi(Acquisition2_High)[[1]], digits = 2),
round(hdi(Acquisition1_Low)[[1]], digits = 2),
round(hdi(Acquisition2_Low)[[1]], digits = 2),
round(hdi(Extinction1_High)[[1]], digits = 2),
round(hdi(Extinction2_High)[[1]], digits = 2),
round(hdi(Extinction1_Low)[[1]], digits = 2),
round(hdi(Extinction2_Low)[[1]], digits = 2))
posterior_means$high_95_CI = c(round(hdi(Baseline1_High)[[2]], digits = 2),
round(hdi(Baseline2_High)[[2]], digits = 2),
round(hdi(Baseline1_Low)[[2]], digits = 2),
round(hdi(Baseline2_Low)[[2]], digits = 2),
round(hdi(Acquisition1_High)[[2]], digits = 2),
round(hdi(Acquisition2_High)[[2]], digits = 2),
round(hdi(Acquisition1_Low)[[2]], digits = 2),
round(hdi(Acquisition2_Low)[[2]], digits = 2),
round(hdi(Extinction1_High)[[2]], digits = 2),
round(hdi(Extinction2_High)[[2]], digits = 2),
round(hdi(Extinction1_Low)[[2]], digits = 2),
round(hdi(Extinction2_Low)[[2]], digits = 2))
posterior_means$low_50_CI = c(round(hdi(Baseline1_High,0.5)[[1]], digits = 2),
round(hdi(Baseline2_High,0.5)[[1]], digits = 2),
round(hdi(Baseline1_Low,0.5)[[1]], digits = 2),
round(hdi(Baseline2_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition1_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition2_High,0.5)[[1]], digits = 2),
round(hdi(Acquisition1_Low,0.5)[[1]], digits = 2),
round(hdi(Acquisition2_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction1_High,0.5)[[1]], digits = 2),
round(hdi(Extinction2_High,0.5)[[1]], digits = 2),
round(hdi(Extinction1_Low,0.5)[[1]], digits = 2),
round(hdi(Extinction2_Low,0.5)[[1]], digits = 2))
posterior_means$high_50_CI = c(round(hdi(Baseline1_High,0.5)[[2]], digits = 2),
round(hdi(Baseline2_High,0.5)[[2]], digits = 2),
round(hdi(Baseline1_Low,0.5)[[2]], digits = 2),
round(hdi(Baseline2_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition1_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition2_High,0.5)[[2]], digits = 2),
round(hdi(Acquisition1_Low,0.5)[[2]], digits = 2),
round(hdi(Acquisition2_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction1_High,0.5)[[2]], digits = 2),
round(hdi(Extinction2_High,0.5)[[2]], digits = 2),
round(hdi(Extinction1_Low,0.5)[[2]], digits = 2),
round(hdi(Extinction2_Low,0.5)[[2]], digits = 2))
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
"High Reward" = "High reward",
"Low Reward" = "Low reward")
posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
"Acquisition1" = "Training1",
"Acquisition2" = "Training2",
"Baseline1" = "Baseline1",
"Baseline1" = "Baseline2",
"Extinction1" = "Test1",
"Extinction2" = "Test2")
# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`Reaction time` = c(555,555,555,555,555,555)
# Raw data
# Prepare the dataset
data.plot = data.final
# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"
# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
"Acq1" = "Training1",
"Acq2" = "Training2",
"Bsln1" = "Baseline1",
"Bsln2" = "Baseline2",
"Ext1" = "Test1",
"Ext2" = "Test2")
data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
"High_Rew" = "High reward",
"Low_Rew" = "Low reward")
data.plot$`Reaction time` = data.plot$Hits.RTs
# Violin + the model
# tiff("RTs_split.tiff", units="in", width=16, height=5, res=800)
ggplot(data=data.plot, aes(x=`Reward Phase`, y=`Reaction time`, label=`Reward Probability`)) +
# violin of the raw data
geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
# individual participants for the raw data
# geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") +
geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
# mean of the model
stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Reaction time`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+
# 95 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_95_CI, ymax=`Reaction time`-`Reaction time`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
# 50 credible intervals of the model
geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_50_CI, ymax=`Reaction time`-`Reaction time`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
scale_x_discrete(limits=c("Baseline1","Baseline2", "Training1","Training2", "Test1","Test2")) +
facet_grid(.~`Reward Probability`) +
scale_y_continuous(breaks=seq(400,750,50)) +
labs(title="Reaction time", y = "Reaction time") +
theme(
axis.line = element_line(size=1, colour = "black"),
panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
panel.border = element_blank(),
panel.background = element_blank(),
plot.title = element_text(size = 12, face = "bold",hjust = 0.5),
text=element_text(colour="black", size = 11),
axis.text.x=element_text(colour="black", size = 11),
axis.text.y=element_text(colour="black", size = 11),
axis.title=element_text(size=12,colour = "black")) +
theme(strip.background =element_rect(fill="white")) +
theme(strip.text = element_text(size = 12))
# dev.off()Table of means across conditions
# Make a table with conditions
posterior_means = as.data.frame(c("Baseline1 High Reward",
"Baseline1 Low Reward",
"Baseline2 High Reward",
"Baseline2 Low Reward",
"Acquisition1 High Reward",
"Acquisition1 Low Reward",
"Acquisition2 High Reward",
"Acquisition2 Low Reward",
"Extinction1 High Reward",
"Extinction1 Low Reward",
"Extinction2 High Reward",
"Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"
posterior_means$Mean = c(paste(round(mean(Baseline1_High), digits = 2), " [", round(hdi(Baseline1_High)[[1]], digits = 2), " ", round(hdi(Baseline1_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline1_Low), digits = 2), " [", round(hdi(Baseline1_Low)[[1]], digits = 2), " ", round(hdi(Baseline1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Baseline2_High), digits = 2), " [", round(hdi(Baseline2_High)[[1]], digits = 2), " ", round(hdi(Baseline2_High)[[2]], digits = 2), "]"),
paste(round(mean(Baseline2_Low), digits = 2), " [", round(hdi(Baseline2_Low)[[1]], digits = 2), " ", round(hdi(Baseline2_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition1_High), digits = 2), " [", round(hdi(Acquisition1_High)[[1]], digits = 2), " ", round(hdi(Acquisition1_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition1_Low), digits = 2), " [", round(hdi(Acquisition1_Low)[[1]], digits = 2), " ", round(hdi(Acquisition1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition2_High), digits = 2), " [", round(hdi(Acquisition2_High)[[1]], digits = 2), " ", round(hdi(Acquisition2_High)[[2]], digits = 2), "]"),
paste(round(mean(Acquisition2_Low), digits = 2), " [", round(hdi(Acquisition2_Low)[[1]], digits = 2), " ", round(hdi(Acquisition2_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction1_High), digits = 2), " [", round(hdi(Extinction1_High)[[1]], digits = 2), " ", round(hdi(Extinction1_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction1_Low), digits = 2), " [", round(hdi(Extinction1_Low)[[1]], digits = 2), " ", round(hdi(Extinction1_Low)[[2]], digits = 2), "]"),
paste(round(mean(Extinction2_High), digits = 2), " [", round(hdi(Extinction2_High)[[1]], digits = 2), " ", round(hdi(Extinction2_High)[[2]], digits = 2), "]"),
paste(round(mean(Extinction2_Low), digits = 2), " [", round(hdi(Extinction2_Low)[[1]], digits = 2), " ", round(hdi(Extinction2_Low)[[2]], digits = 2), "]"))
names(posterior_means)[2] = "Mean [HDI]"
posterior_means = posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
kable(posterior_means, caption = "Means per condition")| Reward Phase | Reward Probability | Mean [HDI] |
|---|---|---|
| Baseline1 | High Reward | 548.63 [ 535.91 560.72 ] |
| Baseline1 | Low Reward | 551.89 [ 538.38 564.15 ] |
| Baseline2 | High Reward | 545.5 [ 532.36 559.34 ] |
| Baseline2 | Low Reward | 556.12 [ 542.01 569.51 ] |
| Acquisition1 | High Reward | 521.46 [ 509.29 534.18 ] |
| Acquisition1 | Low Reward | 542.74 [ 529.9 555.54 ] |
| Acquisition2 | High Reward | 528.94 [ 516.56 541.59 ] |
| Acquisition2 | Low Reward | 534.7 [ 521.78 547.53 ] |
| Extinction1 | High Reward | 528.27 [ 514.43 542.01 ] |
| Extinction1 | Low Reward | 536.87 [ 522 552.32 ] |
| Extinction2 | High Reward | 528.92 [ 515.37 542.04 ] |
| Extinction2 | Low Reward | 542.76 [ 528.68 557.15 ] |
Inference about the best model
Check the difference between baseline1 vs. baseline2 in high reward
Diff_Bsln1_Bsln2_High = Baseline2_High - Baseline1_High
plotPost(Diff_Bsln1_Bsln2_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline1 vs. baseline2 in low reward
Diff_Bsln1_Bsln2_Low = Baseline2_Low - Baseline1_Low
plotPost(Diff_Bsln1_Bsln2_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline2 vs. acquisition1 in high reward
Diff_Bsln2_Acq1_High = Baseline2_High - Acquisition1_High
plotPost(Diff_Bsln2_Acq1_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)
Check the difference between baseline2 vs. acquisition1 in low reward
Diff_Bsln2_Acq1_Low = Acquisition1_Low - Baseline2_Low
plotPost(Diff_Bsln2_Acq1_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)